The aim of this notebook is to compare the annotations between slan+Mono/macro/DC vs slan+ cells.
library(Seurat)
library(tidyverse)
library(harmony)
library(here)
library(glue)
library(pheatmap2)
# Colors
color_palette <- c("#1CFFCE", "#90AD1C", "#C075A6", "#85660D", "#5A5156", "#AA0DFE",
"#F8A19F", "#F7E1A0", "#1C8356", "#FEAF16", "#822E1C", "#C4451C",
"#1CBE4F", "#325A9B", "#F6222E", "#FE00FA", "#FBE426", "#16FF32",
"black", "#3283FE", "#B00068", "#DEA0FD", "#B10DA1", "#E4E1E3",
"#90AD1C", "#FE00FA", "#85660D", "#3B00FB", "#822E1C", "coral2",
"#1CFFCE", "#1CBE4F", "#3283FE", "#FBE426", "#F7E1A0", "#325A9B",
"#2ED9FF", "#B5EFB5", "#5A5156", "#DEA0FD", "#FEAF16", "#683B79",
"#B10DA1", "#1C7F93", "#F8A19F", "dark orange", "#FEAF16", "#FBE426",
"Brown")
# Functions
process_seurat <- function(seurat_obj, group_by_vars = "gem_id") {
seurat_obj <- seurat_obj %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
RunHarmony(group.by.vars = group_by_vars) %>%
RunUMAP(reduction = "harmony", dims = 1:20)
seurat_obj
}
slan_cells <- readRDS(here("scRNA-seq/results/R_objects/slancytes_revision/3.1.seurat_slan_cells_tonsil_classified_deep_annotation_experiments_2_and_3.rds"))
unfiltered <- readRDS(here("scRNA-seq/results/R_objects/slancytes_revision/1.3.seurat_slan_cells_tonsil_filtered_experiments_2_and_3.rds"))
table(slan_cells$gem_id)
##
## q0qp6qyu_nkrug5p0 rssxkrxg_edqhwe1o
## 914 741
slan_cells <- process_seurat(slan_cells)
slan_cells$library_name2 <- str_remove(slan_cells$library_name, "BCLL_10_")
DimPlot(slan_cells, group.by = "library_name2") +
theme(
legend.text = element_text(size = 15),
legend.position = "top",
axis.text = element_blank(),
axis.ticks = element_blank(),
plot.title = element_blank()
)
slan_markers <- read.delim(here("data/raw_data_figures/BA2020Faseb_SLANmarkers.txt"))
sel_genes <- slan_markers$Gene[slan_markers$Gene %in% rownames(slan_cells)]
slan_cells <- AddModuleScore(
slan_cells,
features = list(
SLAN = slan_markers$Gene[slan_markers$CellType == "SLAN"],
DC = slan_markers$Gene[slan_markers$CellType == "DC"],
MAC = slan_markers$Gene[slan_markers$CellType == "MAC"]
),
name = c("SLAN", "DC", "MAC")
)
FeaturePlot(slan_cells, c("SLAN1", "DC2", "MAC3")) &
scale_color_viridis_c(option = "magma") &
theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
)
DimPlot(slan_cells, group.by = "annotation_20220619", cols = color_palette)
FeaturePlot(slan_cells, "annotation_probability") &
scale_color_viridis_c(option = "magma")
DimPlot(slan_cells, group.by = "annotation_20220619", cols = color_palette)
slan_cells@meta.data %>%
group_by(annotation_20220619, library_name) %>%
summarize(n_cells = n()) %>%
group_by(library_name) %>%
mutate(total_cells = sum(n_cells), pct_cells = n_cells / total_cells * 100) %>%
ggplot(aes(library_name, pct_cells, fill = annotation_20220619)) +
geom_col() +
labs(y = "Percentage of cells (%)") +
scale_fill_manual(values = color_palette) +
theme_classic() +
theme(axis.title.x = element_blank(),
legend.title = element_blank())
VlnPlot(slan_cells, c("SLAN1", "DC2", "MAC3"), group.by = "library_name")
Remove lingering poor-quality cells/doublets
slan_cells$doublet_score <- unfiltered$doublet_score[colnames(slan_cells)]
FeaturePlot(slan_cells, features = "doublet_score") &
scale_color_viridis_c(option = "magma")
slan_cells <- FindNeighbors(slan_cells, reduction = "harmony", dims = 1:20, k.param = 15)
slan_cells <- FindClusters(slan_cells, resolution = 1)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1655
## Number of edges: 51500
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7497
## Number of communities: 14
## Elapsed time: 0 seconds
markers_9 <- FindMarkers(slan_cells, ident.1 = "9", ident.2 = "1", only.pos = TRUE, logfc.threshold = 0.75)
DT::datatable(markers_9, options = list(scrollX = TRUE))
markers_3 <- FindMarkers(slan_cells, ident.1 = "3", only.pos = TRUE, logfc.threshold = 0.75)
DT::datatable(markers_3, options = list(scrollX = TRUE))
DimPlot(slan_cells)
VlnPlot(slan_cells, "nFeature_RNA", pt.size = 0) & NoLegend()
slan_cells <- slan_cells[, slan_cells$seurat_clusters != "9"]
DimPlot(slan_cells)
slan_cells <- process_seurat(slan_cells)
DimPlot(slan_cells)
slan_cells <- FindNeighbors(slan_cells, reduction = "harmony", dims = 1:20)
slan_cells <- FindClusters(slan_cells, resolution = 0.15)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 1571
## Number of edges: 54794
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9138
## Number of communities: 5
## Elapsed time: 0 seconds
DimPlot(slan_cells)
slan_sub <- slan_cells[, slan_cells$seurat_clusters == "1" & slan_cells$library_name2 == "SLANCYTES"]
slan_sub <- slan_sub %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA() %>%
# RunHarmony(group.by.vars = "gem_id") %>%
RunUMAP(reduction = "pca", dims = 1:20)
DimPlot(slan_sub, group.by = "gem_id")
FeaturePlot(slan_sub, c("SELENOP", "C1QA", "MMP9", "MAF"))
# Heatmap
goi_slan <- c("SELENOP", "APOC1", "APOE", "FCGR2A", "MMP9", "MMP12",
"C1QA", "C1QB", "HLA-DRA", "HLA-DRB1")
slan_sub <- FindNeighbors(slan_sub, reduction = "pca", dims = 1:20, k.param = 7)
slan_sub <- FindClusters(slan_sub, resolution = 1)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 370
## Number of edges: 8728
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.5370
## Number of communities: 6
## Elapsed time: 0 seconds
DimPlot(slan_sub)
slan_sub$subcluster <- case_when(
slan_sub$seurat_clusters %in% c("0", "1") ~ "0",
slan_sub$seurat_clusters %in% c("3", "4") ~ "1",
slan_sub$seurat_clusters == "2" ~ "2",
slan_sub$seurat_clusters == "5" ~ "3"
)
barcodes <- slan_sub@meta.data %>%
rownames_to_column("barcode") %>%
arrange(subcluster) %>%
pull("barcode")
input_mat <- as.matrix(slan_sub[["RNA"]]@data[goi_slan, barcodes])
pheatmap2(
input_mat,
scale = "none",
show_rownames = TRUE,
show_colnames = FALSE,
border_color = NA,
cluster_rows = FALSE,
cluster_cols = FALSE,
fontsize_row = 6
)
saveRDS(slan_cells, here("scRNA-seq/results/R_objects/slancytes_revision/4.1.seurat_slan_cells_tonsil_dimensionality_reduction_experiments_2_and_3.rds"))
saveRDS(slan_sub, here("scRNA-seq/results/R_objects/slancytes_revision/4.2.seurat_slan_cells_tonsil_subset_slancytes_experiments_2_and_3.rds"))
sessionInfo()
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=es_ES.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] pheatmap2_0.1.0 fastcluster_1.2.3 parallelDist_0.2.6 gtable_0.3.1 scales_1.2.1 RColorBrewer_1.1-3 glue_1.6.2 here_1.0.1 harmony_0.1.1 Rcpp_1.0.10 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.0 purrr_1.0.1 readr_2.1.3 tidyr_1.3.0 tibble_3.1.8 ggplot2_3.4.0 tidyverse_1.3.2 SeuratObject_4.1.3 Seurat_4.3.0 BiocStyle_2.26.0
##
## loaded via a namespace (and not attached):
## [1] readxl_1.4.1 backports_1.4.1 plyr_1.8.8 igraph_1.3.5 lazyeval_0.2.2 sp_1.6-0 splines_4.2.0 crosstalk_1.2.0 listenv_0.9.0 scattermore_0.8 digest_0.6.31 htmltools_0.5.4 fansi_1.0.4 magrittr_2.0.3 tensor_1.5 googlesheets4_1.0.1 cluster_2.1.4 ROCR_1.0-11 limma_3.54.1 tzdb_0.3.0 globals_0.16.2 modelr_0.1.10 RcppParallel_5.1.6 matrixStats_0.63.0 timechange_0.2.0 spatstat.sparse_3.0-0 colorspace_2.1-0 rvest_1.0.3 ggrepel_0.9.3 haven_2.5.1 xfun_0.37 crayon_1.5.2 jsonlite_1.8.4 progressr_0.13.0 spatstat.data_3.0-0 survival_3.5-0 zoo_1.8-11 polyclip_1.10-4 gargle_1.3.0 leiden_0.4.3 future.apply_1.10.0 abind_1.4-5 DBI_1.1.3 spatstat.random_3.1-3 miniUI_0.1.1.1 viridisLite_0.4.1 xtable_1.8-4 reticulate_1.28 DT_0.27 htmlwidgets_1.6.1 httr_1.4.4
## [52] ellipsis_0.3.2 ica_1.0-3 farver_2.1.1 pkgconfig_2.0.3 sass_0.4.5 uwot_0.1.14 dbplyr_2.3.0 deldir_1.0-6 utf8_1.2.3 labeling_0.4.2 tidyselect_1.2.0 rlang_1.0.6 reshape2_1.4.4 later_1.3.0 munsell_0.5.0 cellranger_1.1.0 tools_4.2.0 cachem_1.0.6 cli_3.6.0 generics_0.1.3 broom_1.0.3 ggridges_0.5.4 evaluate_0.20 fastmap_1.1.0 yaml_2.3.7 goftest_1.2-3 knitr_1.42 fs_1.6.0 fitdistrplus_1.1-8 RANN_2.6.1 pbapply_1.7-0 future_1.31.0 nlme_3.1-162 mime_0.12 ggrastr_1.0.1 xml2_1.3.3 compiler_4.2.0 rstudioapi_0.14 beeswarm_0.4.0 plotly_4.10.1 png_0.1-8 spatstat.utils_3.0-1 reprex_2.0.2 bslib_0.4.2 stringi_1.7.12 highr_0.10 lattice_0.20-45 Matrix_1.5-3 vctrs_0.5.2 pillar_1.8.1 lifecycle_1.0.3
## [103] BiocManager_1.30.19 spatstat.geom_3.0-6 lmtest_0.9-40 jquerylib_0.1.4 RcppAnnoy_0.0.20 data.table_1.14.6 cowplot_1.1.1 irlba_2.3.5.1 httpuv_1.6.8 patchwork_1.1.2 R6_2.5.1 bookdown_0.32 promises_1.2.0.1 KernSmooth_2.23-20 gridExtra_2.3 vipor_0.4.5 parallelly_1.34.0 codetools_0.2-19 MASS_7.3-58.2 assertthat_0.2.1 rprojroot_2.0.3 withr_2.5.0 sctransform_0.3.5 parallel_4.2.0 hms_1.1.2 rmarkdown_2.20 googledrive_2.0.0 Rtsne_0.16 spatstat.explore_3.0-6 shiny_1.7.4 lubridate_1.9.1 ggbeeswarm_0.7.1